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O . The deepest space and ground-based observations find metal-enriched galax- 

ies at cosmic times when the Universe was < 1 Gyr old. These stellar popu- 
lations had to be preceded by the metal-free first stars, Population III. Recent 
^ ■ cosmic microwave background polarization measurements indicate that stars 

started forming early when the Universe was 1 200 Myr old. Theoretically it is 
'. now thought that Population III stars were significantly more massive than the 

^ I present metal-rich stellar populations. While such sources will not be individ- 

^ ; ually detectable by existing or planned telescopes, they would have produced 

^ ■ significant cosmic infrared background radiation in the near-infrared, whose 

fluctuations reflect the conditions in the primordial density field. Here we 
report a measurement of diffuse flux fluctuations after removing foreground 
i stars and galaxies. The anisotropics exceed the instrument noise and the more 

c3 ; local foregrounds and can be attributed to emission from massive Population 

III stars, providing observational evidence of an era dominated by these ob- 
jects. 

The cosmic infrared background (CIB) is generated by emission from luminous objects dur- 
ing the entire history of the Universe including epochs at which discrete objects are inaccessible 
to current telescopic studies [Dill. With new powerful telescopes individual galaxies are now 
found out to redshifts of 2 > 5 — 7, but the period preceding that of the galaxies seen in the Hubble 
Ultra-Deep Field (UDF) remains largely unexplored. The UDF galaxies and even the highest 
z quasars (z^6.5) appear to consist of "ordinary" metal-enriched Population I and II stars, 
suggesting that the metal-free stars of Population III lived at still earlier epochs. Large-scale 
polarization of the cosmic microwave background (CMB) |3J suggests that first stars formed 
early at 2; ~ 20. 
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Current theory predicts that Population III objects were very massive stars with mass > 
IOOM0 0161133. They should have produced a significant diffuse background iQ with an 
intensity almost independent of the details of their mass function |13|. Because much of the 
emission longward of the rest-frame Lyman limit from these epochs is now shifted into the 
near-IR (NIR), these stars could be responsible for producing much, or all, of the observed NIR 
CIB excess over that from normal galaxy populations (see detailed review in Q and refs. lEl 
|9l[T0l[TTl[T2l[T3l|). Two groups llT2l[r3ll have recently suggested that this emission should have 
a distinct angular spectrum of anisotropics, which can be measured if the contributions from 
the ordinary (metal-rich) galaxies and foregrounds can be isolated and removed, and provide an 
indication of the era made predominantly of the massive Populations III stars. 

Measuring CIB anisotropics from objects at high z is difficult because the spatial fluctua- 
tions are small and can be hidden by the contributions of ordinary low-^ galaxies as well as 
instrument noise and systematic errors in the data. Previous attempts to measure the structure 
of the CIB in the near-IR on the degree scale with C05£/DIRBE |21 1, /i?r5/NIRS |29| were 
limited in sensitivity because of the remaining contributions from brighter galaxies in the large 
beams. Analysis of 2MASS data at 1.25, 1.65 and 2.2 yum with ~ 2 arcsec resolution ll33ll34ll 
allowed removal of foreground galaxies to a K magnitude of 19-19.5 (m^^ ~ 21 or 15 /xJy) 
and led to measurements of CIB fluctuations at 1.25 to 2.2 yum at sub-arcminute angular scales. 
These studies reported fluctuations in excess of that expected from the observed galaxy popu- 
lations, although their accurate interpretation in terms of high-2; contributions is difficult due to 
foreground galaxies and non-optimal angular scales. 

We used data from deep exposure data obtained with 5p«Yzer/IRAC lfT5lfT9ll in four channels 
(Channels 1-4 correspond to wavelengths of 3.6, 4.5, 5.8 and 8 /um respectively) in attempting to 
uncover this signal. We find significant CIB anisotropics after subtracting galaxies substantially 
fainter than was possible in prior studies, i.e. down to mAB ~ 22 — 25. The angular power 
spectrum of the anisotropics is significantly different from that expected from Solar System and 
Galactic sources, its amplitude is much larger than what is expected from the remaining faint 
galaxies, and can reasonably be attributed to the diffuse light from the Population III era. 

Assembly and reduction of data sets. The primary data set used here is the deep IRAC 
observation of a ~ 12' x 6' region around QSO HS 1700+6416 obtained by the Spitzer instru- 
ment team lfT?l [T51 . In addition, the data from two auxiliary fields with shallower exposures 
(HZF and EGS) were analyzed. Relevant data characterising the fields are listed in Table 1 . For 
this analysis the raw data were reduced using a least-squares self-calibration method lITSl : the 
processing is described in more detail in the Supplementary Information (hereafter SI). 

The random noise level of the maps was computed from two subsets (A, B) containing 
the odd and even numbered frames, respectively, of the observing sequence. The A and B 
subsets were observed nearly simultaneously and with similar dither patterns and exposures. 
The difference between maps generated from the A and B subsets should eliminate true celestial 
sources and stable instrumental effects and reflect only the random noise of the observations. 

Analysis of the background fluctuations must be preceded by steps that eliminate the fore- 
ground Galactic stars and the galaxies bright enough to be individually resolved. The primary 
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Table 1 : Analyzed fields, their coordinates, exposure times, limiting magnitudes, and pixel 
scales 



Region 


QSO 1700 (Ch 1-4) 


HZF(Ch 1-3) 


HZF (Ch 4) 


EGS (Ch 1-4) 


{a, 6) 


(255.3, 64.2) 


(136.0, 11.6) 


(285.7, -17.6) 


(215.5,53.3) 


(l, &)Gal 


(94.4, 36.1) 


(217.5, 34.6) 


(18.4, -10.4) 


(96.5, 58.9) 


(A,/5)ec1 


(194.3, 83.5) 


(135.0, -4.9) 


(285.0, 5.0) 


(179.9, 60.9) 


{tohs){hrs) 


7.8, 7.8, 7.8, 9.2 


0.5, 0.5, 0.5 


0.7 


1.4, 1.4, 1.4, 1.4 


^Vcga,lim 


22.5, 20.5, 18.25, 17.5 


21.5, 19.5, 17 


14.5 


22.5, 20.75, 18.5, 17.75 


Pixel scale (") 


0.6 


1.2 


1.2 


1.2 


Field size (pix) 


1, 152 X 512 


576 X 256 


576 X 256 


640 X 384 



QSO 1700 was observed during In-Orbit Checkout (IOC) with eight AORs (Astronomical Ob- 
servation Requests) which used various dither patterns and 200-sec frame times, except for two 
which used 100-sec frame times (AOR ID numbers = 7127552, 7127808, 7128064, 7128320, 
7128576, 7475968, 7476224, 7476480). Because of the focal plane offset between the shared 
3.6/5.8 iim and 4.5/8 iim fields of view, each of these pairs observes separate fields which 
have a common overlap of ~ 5' x 5' at all four wavelengths. To provide contrast for the self- 
calibration algorithm, data from a high-zodiacal light brightness field (HZF) was co-processed 
for each channel. For 3.6, 4.5, and 5.8 fim, the nearest suitable (200-second frame time) data 
were observed later during IOC (AORID number = 8080896). For the 8 /xm data, because 
nominal 100 and 200-second frame times are split into pairs and quartets of 50-second frames, 
more nearly contemporaneous observations earlier during IOC were used (AORID number = 
6849280). For this work we only self-calibrated the six 200-second AORs for 3.6 - 5.8 /im, but 
at 8 fim all eight AORs were used since all produce data with 50-second frame times. Observa- 
tions of the extended Groth Strip area (EGS) provide an additional deep data set for verification 
of our results. Separate exposure times and limiting magnitudes apply for Channels 1-4 respec- 
tively. The maps for the main QSO 1700 field are shown in SI (the supplementary information 
available at Nature online). 
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means of removing these sources is an iterative clipping algorithm which zeros all pixels (and 
a fixed number of neighboring pixels, A^mask x ^mask) which are more than a chosen factor, 
A^'cut, above the la RMS variation in the clipped surface brightness. This must be restricted to 
relatively high Ncut to (a) leave enough area for a robust Fourier analysis of the map, and (b) 
avoid clipping into the background fluctuation distribution. This means that faint sources, the 
faint outer portions of resolved galaxies, and the faint wings of the point source response func- 
tion around bright stars cannot be clipped adequately. To remove these low surface brightness 
sources we used a CLEAN algorithm [ 16 1 to model the entire field in each channel. This model, 
convolved with the full IRAC point spread function (PSF), was subtracted from the undipped 
regions of the map. SI provides details on this process and illustrations of the clipped and 
model-subtracted images. The final step is the fitting and removal of the 0th and 1st order com- 
ponents of the background in the undipped regions. This is done to minimize power spectrum 
artifacts due to the clipping, and because the lack of an absolute flux reference measurement 
and observing constraints prohibit unambiguous determination of these components. With this 
subtraction, the images represent the fluctuation fields, 5F{x), rather than the absolute intensity. 
For each observed field these steps were carried out for images derived from the full data set 
and from the A and B subsets. 

Power spectrum computation and analysis. For each channel, we calculate the power 
spectra of the fluctuations as a quantitative means of characterizing their scale and amplitude. 
The power that remains after subtraction of the random noise component is insensitive to the de- 
tails of the source clipping, and is statistically correlated between channels. This confirms that 
the fluctuation signal does have a celestial origin. The shape and amplitude of the power spec- 
tra are not consistent with significant contributions from the cirrus of the interstellar medium 
(except perhaps at 8/im) or from the zodiacal light from local interplanetary dust. 

The fluctuation field, 5F{x), was weighted by the observation time in each pixel, w{x) oc 
tobs(^c), and its Fourier transform, f{q) = J 5F{x)w{x) ex\){—ix-q)(P'x calculated using the 
fast Fourier transform. (The weighting is necessary to minimize the noise variations across the 
image, but we also performed the same analysis without it and verified that the weight adds no 
structure to the resultant power spectrum. This is because the weights are relatively flat across 
the image). The power spectrum is ^2(9) = with the average taken over all the 

Fourier elements Nq corresponding to the given q. A typical flux fluctuation is \J P2{q) / 27: 
on the angular scale of wavelength 2n / q. In SI we show the final power spectrum of the diffuse 
flux fluctuations, Ps{q), and the noise, -PAr(q'), of the datasets. We find significant excess of the 
large-scale fluctuations over the instrument noise. 

Fig. [U shows the excess power spectrum, Psiq) — PNiQ), of the diffuse light after the in- 
strument noise has been subtracted. There is a clear positive residual whose power spectrum is 
significantly different from white noise and has substantial correlations all the way to the largest 
scales probed. Possible sources of these large-scale correlations can be artifacts of the analysis 
procedure (e.g. clipping and FFT), instrumental artifacts, local Solar System or Galactic emis- 
sion, or relatively nearby extragalactic sources and/or more distant cosmological sources. In 
what follows we discuss their relative contributions. 
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Table 2: Cross-correlation of fluctuations in units of that of random sample IZ/gti for 
various clipping thresholds 



Channels 


iVcut = 4 


A^cut = 2 


1 : 2 


52 


12 


1 : 3 


7 


0.6 


1 : 4 


10 
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These correlations were evaluated for an evenly covered region of 512^ pixels common to all 4 
channels for the QSO 1700 field. For random uncorrelated samples, 51,^2, of A^pixeis the cor- 
relation coefficient, IZ = {^1^2) /[{^D^^DY^"^^ should be zero with dispersion a-ji = N'^^^"^^. 
Whereas, for A^mask = 3, we find that down to even Ncut = 2, when only 6% of the pixels 
remain, the correlations between the channels remain statistically significant. Simple simula- 
tions containing the appropriate levels of the instrument noise and a power law component of 
the CIB gave somewhat larger mean values of TZ with the measured values lying within a 95% 
confidence level; incorporating the possibility that not all of the remaining ordinary galaxies are 
the same at each wavelength would however reduce the mean values of TZ. 



Residual emission from the wings of incompletely clipped sources can give rise to spurious 
fluctuations, but the power spectrum of these fluctuations should depend on the clipping param- 
eters iVcut and A^mask- We tested the contributions from these residuals in various ways. For a 
given Ncut we varied A'mask from 3 to 7 significantly reducing any residual wings and increasing 
the fraction of the clipped pixels, but found negligible a few percent) variations in the final 
P2{q)- We also clipped down to progressively lower values of A^cut- For Ncut ^ 3.5 too few pixels 
remain for robust Fourier analysis, so in these cases we computed C (9), the correlation function 
of the diffuse emission, related to the power spectrum by a 1-D Legendre transformation. It is 
consistent with the fluctuations in Fig. [T]as discussed and shown in SI. 

Instrument noise contributions to the fluctuations were evaluated as the power spectra of 
^{5Fa — SFb) using the A — B subset images. As shown in SI, random instrument noise has 
an approximately white spectrum. The power spectra of the final datasets have a much larger 
amplitude than the noise (especially in Channels 1 and 2) until the convolution with the beam 
tapers off the signals from the sky at the smallest angular scales. The instrument noise spectra 
are unaffected by the beam and are uncorrelated from channel to channel, as expected. However, 
the full data set fluctuation fields show statistically significant correlations between channels as 
shown in Table 2. This means that we see the same fluctuation field in addition to (different) 
noise in all four channels. SI describes tests to assess the contributions of possible instrumental 
systematic errors. They indicate that systematic instrumental effects are unlikely to lead to the 
signal shown in Fig. [H 

The best assessment of zodiacal light contributions to the power spectrum comes from the 
examination of EGS observations taken at two epochs 6 months apart. Because any anisotropics 
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in the zodiacal light cloud will not remain fixed in celestial coordinates over this interval, the 
difference in the fluctuation fields at these two epochs should eliminate Galactic and extragalac- 
tic signals and yield a power spectrum of the zodiacal light fluctuations added to the instrument 
noise and possible systematic errors. The fluctuation levels of these difference maps set an up- 
per limit on the zodiacal light contribution of < 0.1 nW/m^/sr at 8 yum. Scaling this result to 
the other channels by interpolating the observed zodiacal light spectrum ll36l leads to zodiacal 
light contributions to the fluctuations that are comfortably below the detections in Fig. [l]in the 
other channels as well. 

Our assessment of the contribution of the infrared cirrus (i.e. interstellar clouds of neutral 
gas and dust) to the power spectra is derived from the 8 /im HZF field which lies at low fecai 
and is visibly contaminated by cirrus. Assuming that the large scale fluctuations in this field 
are due to the cirrus, relative fluctuations of the 8 fxm cirrus are ~ 1% of the mean cirrus 
flux level. A similar level of relative cirrus fluctuations in the QSO 1700 field would have an 
amplitude of ~ 0.3 nW/m^/sr. This is not significantly lower than the 8 fim amplitude observed 
in Fig. [H Therefore, we cannot presently eliminate the possibility that the fluctuations at 8 /im 
are dominated by cirrus. However, the spectrum of the ISM emission should drop sharply at 
shorter wavelengths as the PAH emission bands that dominate at 8 fim become less significant. 
Given this estimate of the cirrus contribution at 8 fim we estimate the amplitude of large-scale 
fluctuations for channels 1 - 3 as ~ 0.03, 0.03, and 0.08 nW/m^/sr, which are well below the 
observed fluctuations in these channels. 

Colour is another important criterion for testing the origins of the signal in Fig. [B The 
cross-correlations between the channels for the common area of the QSO 1700 field are statis- 
tically significant and they also strengthen at the larger scales where the noise contribution is 
smaller, as verified by smoothing the maps. These significant cross-correlations allow us to ex- 
amine the corresponding colours of the fluctuations. We made several estimates of the colours, 
Pni, between channels n and 1: a) as the square root of ratio of the power spectra, and b) as 
Pni{d) = {5n{x)5i{x + 9)) / {5i{x)Si{x + 9)) evaluated over the 512^ pixel field common to 
all four channels, where a further consistency check comes from comparison between P21 and 
P4:i/Pi2- For channels 1 and 2 these estimates are shown in Fig. Eland appear roughly inde- 
pendent of angular scale and mutually consistent. The instrument noise is too large to enable 
robust colour estimates involving channel 3 (5.8 fim). These derived colours indicate that the 
fluctuation signal in Fig. [l]has an energy distribution that is approximately flat to slowly rising 
with wavelength in ul^. The energy spectrum rules out contributions from remaining Galactic 
stars, but probably cannot be used to distinguish between ordinary galaxies and Pop III objects 
without additional detailed modeling of both sets of sources. 

CIB fluctuations from extragalactic sources. There are two extragalactic classes of con- 
tributors to CIB fluctuations: "ordinary" galaxies containing normal stellar Populations I and 
II, and the objects that preceded them. Population III stars. The square of the CIB fluctuation 
in band v produced by cosmological sources that existed over time period At is given by the 
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power- spectrum version of the Limber equation (see |l2| and SI): 




2 



A2 



[qdj^] z)dt 



(1) 



where v' = z/(l + z) is the rest frequency of the emitters, (Ia is the comoving angular diameter 
distance and 



is the fluctuation in the number of sources within a volume k ^cAt. The fluctuation of the CIB 



the suitably averaged effective redshift. The relative CIB fluctuation is a suitably averaged fluc- 
tuation in the source counts over a cylinder of radius q^^dA{{z)) and length cAt. Three things 
would lead to larger fluctuations: 1) a population that, after removing constituents brighter than 
some limit, leaves a substantial mean CIB flux (increase Fcm), 2) populations that existed for a 
shorter time (increase (A) by decreasing At), and 3) populations that formed out of rare peaks 
of the underlying density field leading to biased and significantly amplified, [27| clustering 
properties (increase (A) by increasing Pi^k)). 

CIB fluctuations from remaining ordinary galaxies. The fluctuations produced by or- 
dinary galaxies contain two components: 1) shot noise from discrete galaxies and 2) galaxy 
clustering from the primeval density field. The amplitude of the first component can be es- 
timated directly from galaxy counts. In order to estimate the contribution from the second 
component we proceed as follows: from galaxy count data we estimate the total CIB flux pro- 
duced by the remaining galaxies fainter than our clipping threshold. Ordinary galaxies occupy 
an era of At^ a few Gyrs. Their present-day clustering pattern on the relevant scales is well 
measured today. The clustering pattern evolution can be extrapolated to earlier times assuming 
the "concordance" ACDM model. These parameters (flux. At, clustering pattern and its evo- 
lution) then allow us to estimate the contribution to the CIB fluctuations via eqs. 11121 Because 
galaxy clustering is weaker at earlier times, an upper limit on the contribution from ordinary 
galaxies is obtained assuming that the clustering at early times remained the same as at 2; = 0. 

The shot noise component contributed by the ordinary galaxies to the CIB angular spectrum 
was estimated directly from galaxy count data: each magnitude bin Am with galaxies of flux 
f{m) would contribute a power spectrum of Psn = /^(m)^^^Am. Fig. [l] shows the shot 
noise component convolved with the IRAC beam for the galaxies at the limiting magnitude in 
Table[T] The limiting magnitudes are qualitative estimates of where the SExtractor |22| number 
counts (in 0.5 magnitude bins) begin to drop sharply due to incompleteness. The shot noise fits 
the observed diffuse light fluctuations well at small angles, but the shot-noise contribution from 
ordinary galaxies cannot make a substantial contribution to the large scale power of the diffuse 
light in Fig. [H 

We assume the "concordance" ACDM model Universe with flat geometry (i^totai=l) dom- 
inated by a cosmological constant, ~ 0.7, with the rest coming from cold dark matter and 




(2) 



on angular scale ~ 27r/g can then be expressed as 6 Fcm 



FciBA{qd/{{z))), where {z) is 
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ordinary baryons in proportions suggested by measurements of the CMB anisotropies H and 
high-2; supernovae ll47l . With P^ik) taken from the ACDM concordance model, the power in 
fluctuations from the clustering of ordinary galaxies depends on the net flux produced by the 
remaining ordinary galaxies (via Fcib) and their typical redshifts (via A(qd^^)). The total flux 
from galaxies fainter than the limiting magnitude in Table 1 was estimated from the Spitzer 
IRAC galaxy counts ( ifTSl ) and is shown in Fig. |3l The flux contributed to the CIB from the 
remaining galaxies is between 0.1 and 0.2 nW/m^/sr and this amplitude is much less than the 
excess CIB at these wavelengths ([2|) indicating that the excess CIB flux is produced by popu- 
lations that are still farther out than the ordinary galaxies remaining in the data. This (low) value 
of the remaining flux can also be derived from the small amplitude of the residual shot-noise 
contribution to the fluctuations in the confusion-limited datasets. Thus, in order to explain the 
amplitudes shown in Fig. [H the remaining ordinary galaxy populations would need to produce 
relative flux fluctuations of order > 100% from their clustering. Without follow-up spectroscopy 
it is difficult to determine with high precision the range of redshifts of the remaining ordinary 
galaxies, but approximate estimates can be made and are summarized in the Figure |3] legend to 

In flat cosmology with il/^ = 0.7 one arcminute subtends a comoving scale between 0.7 
and 1.6 h~^Mpc at z between 1 and 5. The present day 3-D power spectrum of galaxy cluster- 
ing, P3,gai(fc), is described well by the concordance ACDM model (|25J, [i26|) and one expects 
that on arcminute scales the density field was in the linear to quasi-linear regime at the red- 
shifts probed by the remaining galaxies between 3.6 and 8 fim. At smaller scales non-linear 
corrections to evolution were computed following (||44||). The resultant CIB fluctuation from 
the remaining ordinary galaxies, producing -Fcm.og = 0.14 nW/m^/sr , times A((x (At)"^/^) 
is shown in Fig. [T]for (2)=!, 3,5 and At=5 Gyr corresponding to the age of the Universe at 
2 ~ 1. For the ACDM model the relative fluctuations in the CIB on arcminute scales would 
be of order (A) ~ (2 — 10) x 10^^(At/5Gyr)^°-^ from galaxies with {z)=\-5 assuming no 
biasing. Combining this with the above values for the diffuse flux from the remaining ordinary 
galaxies would lead to 5Ft (1 — 2) x 10^^ nW/m^/sr in all the channels. While biasing may 
increase the relative fluctuations, with reasonable bias factors for galaxies lying at 2; ~ (a few), 
the diffuse light fluctuations would still be very small compared to those in Fig. [T]and are un- 
likely to account for fluctuations of amplitude ~ (0.1 — 0.5) nW/m^/sr at arcminute scales. 
An upper limit on the CIB fluctuations can be evaluated assuming the same clustering pattern 
for the remaining galaxies as at the present epoch, z = 0, i.e. that their 2-point correlation func- 
tion is given by ^ = (r/r^,)^^ '' with = 5.5/i~^Mpc f4^; its A times -FciB.og is also shown 
in FigUJand is much below the signal we measure. Conversely, one can evaluate the cluster- 
ing strength needed to account for the observed fluctuations: a 100% relative fluctuation from 
galaxies at ztl clustered with ^ = (r/r*)^^'^ would require ^ 25/i^^Mpc. This corresponds 
to the effective bias factor Z 5 which is significantly higher than expected from gravitational 
clustering evolution in the ACDM Universe [351 . In fact, direct measurements of clustering for 
relatively nearby IRAC galaxies with flux > 32/iJy at 3.6 fim (~5 magnitudes brighter than the 
limit reached for remaining galaxies in the QSO 1700 field) find the projected 2-point cor- 
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relation amplitude on ~arcmin scales of u'>32^jy(6'ii 1') < 4 x 10^^ corresponding to relative 
fluctuations in CIB from these galaxies of amplitude ~ y/wtO.2. Fainter galaxies are expected 
to have an even lower correlation amplitude. 

The QSO 1700 field seems to contain an overdensity of Lyman-break galaxies at z ~ 2.3 
ll45ll . which could lead to a larger A and CIB fluctuations for this region. However, we see 
similar levels of fluctuations for the other fields located at very different parts of the sky making 
it unlikely that the overdensity claimed in ref. ll45l can account for our signal. When account is 
made of the different shot noise levels from the remaining galaxies, the fluctuations seen in the 
different fields have consistent power spectra within the statistical uncertainties. (An exception 
being channel 4 HZF data, located at low fecai, and clearly dominated by Galactic cirrus). 

CIB fluctuations from the Population III era. Population III at z ~ 10-30 are expected to 
precede ordinary galaxy populations. One can expect on fairly general grounds lfT3ll that, if mas- 
sive, they would contribute significantly to the NIR CIB, both its mean level and anisotropics. 
Intuitive reasons are discussed in [2|, but as eqs. 1 1121 shows they are mostly related to: 1) if 
massive. Population III were very efficient light emitters, 2) their era likely lasted a shorter 
time. At, than that of the ordinary galaxies, leading to larger A in eq|2l and 3) they should have 
formed out of high peaks of the density field whose correlation function is strongly amplified. 
The near-IR also probes UV to visible parts of the electromagnetic spectrum at z ~ 10 — 30, 
where most of their emission is produced [1 IJ. 

The amplitude of the CIB anisotropics remaining in the present data implies that the re- 
maining CIB originates from still fainter objects. Can the observed amplitudes of fluctuations 
in Fig. [T]be accounted for by energetic sources at high 2? Population III stars can produce 
significant NIR CIB levels [11|, ^1 nW/m^/sr and e.g. the NIR CIB excess over that from 
ordinary galaxies at 3.6 fim is 8.7 ±3.1 nW/m^/sr [2J. Thus Population III objects would re- 
quire smaller relative CIB fluctuations. Because individual Population III systems are small, yet 
numerous, the shot-noise component of the CIB from them is small and, in any case, is already 
absorbed in the shot-noise shown in Fig. [H Additionally, their (A) would be amplified by the 
much shorter At and (significant) biasing, and Population III contribution would dominate the 
diffuse light fluctuations in Fig. [B Detailed theoretical interpretation of the results in terms 
of the Population III era models exceeds the scope of this article, but qualitative comparison 
can be made by estimating the typical value of the relative CIB fluctuation. A, corresponding 
to that era. The fraction of the Population III haloes was calculated assuming that they form 
from the ACDM density field in haloes where the virial temperature T^iy- > 2000 K to enable 
efficient molecular hydrogen cooling ([38|). Biasing was treated using the gravitational cluster- 
ing prescription from |39| in the ACDM model. (Non-linear evolution effects are small on the 
angular scales and redshifts probed here.) Assuming At=300 Myr, which corresponds to the 
age of the Universe at z=lO, we get typical values of A ~ 0.1-0.2 at z between 10 and 20. This 
order-of-magnitude evaluation shows that the levels of ~0. 1-0.3 nW/m^/sr at 3.6 to 8 fim on 
arcminute scales can be accounted for by Population III emissions if their total flux contribution 
is > 1 nW/m^/sr which is reasonable as discussed earlier. 

Earlier studies of CIB fluctuations II2TI . 113311 . 11291 contained significant contributions from 
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relatively bright galaxies making it difficult to isolate the possible CIB fluctuations from the 
very early times. The contribution to CIB fluctuations from remaining galaxies is a function 
of the limiting magnitude below which galaxies are removed. With the Spitzer IRAC data we 
could identify and remove galaxies to very faint limits of flux ^0.3/iJy. This limit is, at last, 
sufficiently low to push the residual contribution from ordinary galaxies along the line-of-sight 
below the level of the excess signal at larger angular scales. If our interpretation is correct and 
the signal we detect comes from Population III located at much higher z, the amplitude of the 
CIB fluctuations on scales where galaxy shot-noise is negligible should remain the same as 
fainter ordinary galaxies are removed with deeper clipping. This is true as far as we can test 
with this data and would certainly be verifiable with longer exposure data. 

At these z the IRAC bands probe the rest-frame wavelengths between 0.2 and 0.8 //m, where 
the energy spectrum of individual Population III object is dominated by free-free emission and 
is a slowly rising function of wavelength [ 1 1 1; this would be consistent with the spectrum of the 
CIB anisotropics in Fig. [l] Near-IR observations at shorter wavelengths would be particularly 
important in confirming the redshifts where the CIB fluctuations originate as there should be a 
significant drop in the fluctuations power at the rest-frame Lyman limit wavelength. 
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Supplementary Information for "Tracing the first stars witli cosmic infrared back- 
ground fluctuations" 

by A. Kashlinsky, R. G. Arendt, J. Mather, S. H. Moseley 

Tliis material presents technical details to support the discussion in the main 
paper. We discuss here data assembly, removal of individual resolved sources, 
computations of the power spectra and correlation functions from the maps 
cleaned off these sources, and evaluation of the Limber equation in the form 
used in the main paper. We also show maps of the cleaned field and present 
the plots with their power spectrum and that of the instrument noise. 

Data processing 

For this analysis the raw data were reduced using a least-squares self-calibration method 
(ref. [18] in main paper). This procedure models the raw data (D^) of each frame as 

= 5" + FP + F1 (SI-1) 

where the index i counts over all pixels of all frames, S"' is the sky intensity at location a, 
is the detector gain at pixel p, is the detector offset at pixel p, and is an offset that can 
differ for each frame and each of the four detector readouts (combined into the index q). This 
self-calibration procedure has advantages over the standard pipeline calibration of the data in 
that the derived detector gain and offsets match the detector at the time of the observation, rather 
than at the time of the calibration observations. The most noticeable difference (especially at 
3.6 and 8 /xm) is that the different AORs are affected by different patterns of residual images 
left by prior observations, which if not properly removed, can leave artifacts in the final images. 

The combined use of low (QSO 1700) and high brightness (HZF) fields removes the gain- 
offset degeneracy in the self-calibration, although since we are interested in the nearly flat 
uniform background, to first order it is more important that the observing strategy included 
dithering to distinguish between detector and sky variations. Because IRAC's calibration does 
not include zero-flux closed-shutter data, a degeneracy remains in the absolute zero point of 
the solution. The sky intensity can be increased by a constant value at all locations a, and 
the offset F^ can be appropriately decreased without affecting the quality of the fit. A de- 
generacy also remains in first-order gradients (or two-dimensional linear backgrounds). The 
parameter can correlate with the dither position to induce false linear gradients in the other 
parameters without affecting the quality of the fit. However higher-order gradients are not de- 
generate. Because of these degeneracies, linear gradients have been fit and subtracted from the 
self-calibrated mosaicked images. Hence our estimates of the power at scales > 5' should be 
treated as lower limits. 

After the self-calibrations had been determined, corrections for "column pull-down" at 3.6 

and 4.5 /im and banding or "row pull-up" at 8 fim were applied to the calibrated individual 
frames before mosaicking into a final image. Both these electronic effects only occur at bright 
sources, and are strictly aligned with the array coordinates and thus would also be minimized 
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by selective filtering in the Fourier domain (discussed below). The column pull-down algorithm 
applied is the "do_pulldown" procedure which is available on the Spitzer Science Center (SSC) 
website as a contributed software tool for post-BCD (Basic Calibrated Data) processing. The 
row pull-up correction used at 8 /xm adjusts each row of each frame as a function of the total 
flux above the 90th-percentile level of that row. The mosaicking built into the self-calibration is 
an interlacing algorithm, in which each datum (D^) is averaged into only a single sky pixel (5"^). 
This means that the random noise in translates into a completely random noise component 
in 5*", but it also means that for sky maps with half-size pixels (0.6") as we use here for the 
QSO 1700 data, the effective coverage for a sky pixel (a) is reduced by a factor of 4. 
Modeling Resolved Sources 

In the full data set mosaics, the foreground sources were identified and measured using 
SExtractor (ref. [22] in main paper). The resulting source catalog is statistically consistent 
(within the Poissonian uncertainties) with the IRAC source counts (ref. [15] in main paper), 
both in terms of numbers of sources and limiting magnitudes. A model image of the foreground 
sky was then constructed by scaling an IRAC point spread function (PSF) image by the inten- 
sity of each source and adding it to an initially blank map at the appropriate location of each 
source. For one version of this, we used the PSFs that were determined during IOC and are 
available as FITS images on the SSC website. However, these PSFs are only mapped within 
a ~ 24" X 24" region, and thus neglect the extended low-level wings of the PSFs. Thus a 
second model was created in which the PSFs applied consisted of the IOC PSF cores smoothly 
grafted onto the extended wings as measured by observations of a bright (saturated) star. ^ This 
wide grafted PSF avoids saturation at the core while tracing the extent out to the full size of the 
array (~ 300" x 300"). However, even with these wide PSFs, these models cannot account for 
the extended emission around many of the galaxies that are resolved by IRAC. Therefore we 
developed a final model, similar to a CLEAN algorithm (ref. [16] in main paper). To construct 
this model for each channel (1) the maximum pixel intensity was located (and saved to a list of 
"clean components"), (2) the wide PSF was scaled to half of this intensity and subtracted from 
the image, (3) this process was iterated 60,000 times, saving intermediate results every 3,000 it- 
erations. The results were examined to identify the appropriate iteration number where the clean 
components become largely uncorrected with the sources in the original image (i.e. they rep- 
resent randomly distributed peaks without apparent convolution by the PSF). This limit occurs 
much more quickly in the longer wavelength channels, where the relative noise level is higher 
and the number of detected sources is fewer. The "model" is then the designated number of 
clean components convolved with the wide PSF. Since this model can represent extended emis- 
sion with groups of clustered components, it is much better at removing the extended emission 
of the resolved galaxies. However, even this model will miss extended emission that is below 
the pixel-to-pixel noise. Such low surface brightness features must be dealt with as part of the 
analysis and interpretation of the power spectra. 

Clipping individual sources 



^Fomalhaut, AOR ID number = 6066432 
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The original maps were clipped of resolved sources (Galactic stars and galaxies) using an 
iterative method developed for DIRBE data analysis (ref. [21] in main paper). Each iteration 
calculates the standard deviation (a) of the image, and then masks all pixels exceeding A'^cut cr 
along with A^mask x ^mask surrouudiug pixels. The procedure is repeated until no pixels exceed 
iVcut cr. This clipping tends to miss emission in the outer parts of extended sources and the 
distant wings of very bright point sources. To help minimize this effect, the clipping (and 
subsequent analysis) is followed by subtraction of the models described above. We also clip 
the model images themselves with the same algorithm at the same level of Ncut with the mask 
derived from the clipping of the actual data field. The final mask was taken to be the combined 
mask from the two clipped maps, the main data and the model map. Whereas the appearance of 
the residual emission around bright sources has decreased dramatically with this extra step, the 
results were practically unaffected by this extra clipping (< a few percent in large scale power). 

Truly robust results must be independent of the clipping parameters to within the variations 
from the signal arising from the populations of sources corresponding to the different clipping 
thresholds, Ncut- We explored the range for the various clipping parameters: the results are 
largely independent of the mask size and in what follows a value of A'^mask = 3 was adopted for 
displaying the final power spectra, which corresponds to ^ 75 % pixels remaining in the clipped 
field at Ncut = 4. For sufficiently low values of Ncut too few pixels remained in the field for 
robust Fourier analysis of the diffuse light distribution. In order to further verify robustness of 
the results we also evaluated the correlation function at lower values of Ncut when significantly 
fewer pixels remained. This simple clipping algorithm proved quite efficient in identifying the 
sources. Comparison with the SExtractor model catalog shows that in channels 1 and 2 our 
clipping algorithm identified over 95% of the sources and in channel 3 and 4 the overlap with 
the SExtractor catalog was nearly 100%. 

Computing power spectrum from clipped maps 

The blanked pixels were assigned 5F=0, thereby not adding power to the eventual power 
spectrum. Prior to computation of the Fourier transform and the correlation function, we sub- 
tracted remaining linear gradients from the maps. The Fourier transform image was cleansed of 
stripes and muxbleed artifacts in the {u,v) transform space by zeroing the power at frequencies 
corresponding to n = ±nqmux{n= 1,2, 3,...) where 7i/qmux = 4 x 1.2" which is the spacing 
of muxbleed artifacts. Power was also set to zero along the u=0 and f =0 axes to eliminate any 
residual signal from the gradient subtraction. The final power spectrum of the residual diffuse 
light was computed by averaging over narrow concentric rings of radius q. The noise was eval- 
uated from the A — B difference maps using the same clipping mask and weighting method as 
the final co-added data. 

Fig. 1ST 1 1 shows the clipped maps for the QSO 1700 field used in the analysis and Fig. IST3I 
plots the power spectra deduced from the clipped maps and their noise. 

Fig. IST2I shows the maps clipped at Ncut = 2 with A^mask = 3. Lowering A'cut leaves huge 
multi-connected holes in the image. This precludes robustly computing Fourier transform of 
the maps; hence we evaluated the autocorrelation function, C{9) = {5F(x + 6)5F{x)) for low 
A'cut, which is shown in Fig. ISI-4[ The figure shows that at every level of A'cut we find the same 
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large scale correlations with C{9) remaining positive to at least ~ 100". (The slight drop in the 
amplitude with lower iVcut is expected as progressively lower peaks of CIB and the instrument 
noise are clipped). Note that the power spectrum, if of cosmological origin, may depend on 
A'cut if the clipping pushes deep enough to start eliminating the high-2; cosmological sources. 
Systematic instrument effects 

Systematic errors associated with the detector self -calibration, i.e. the gain (or flat field) and 
the offset (or dark frame), should produce a pattern on the sky which is the fixed pattern of the 
error convolved with the dither pattern of the pointings for each AOR. Since we do not know the 
pattern of the systematic error, we created mosaics from the detector flat fields, from the detector 
dark frames, and from frames consisting of a single Gaussian point source at a fixed position 
on the detector. Visually, the sort of structure that appears in these systematic error images 
does not look like the actual sky images. This is confirmed quantitatively by performing the 
power spectrum analysis on these images, which in all cases show fluctuations with a distinctly 
different spectrum from those of the real data. A different test for possible instrumental effects 
is to check the isotropy of the fluctuations. Instrumental effects (such as muxbleed, column 
pull-down, etc.) are likely to be associated with preferred directions on the array and on the 
sky as well since the position angle of the detector was roughly constant for the duration of the 
observations. As a test of isotropy for each channel we calculated the power spectra in four 45°- 
wide azimuthal sectors of the Fourier transform's (ti, v) plane centered along the lines v = 0, 
V — u,u — and v — —u. The sectors along the axes would be expected to be more strongly 
affected by instrumental effects, but we find that there are no large-scale variations between the 
power spectra in the different sectors and only ^5% variations in the overall amplitudes. Stray 
light is a systematic error which remains relatively fixed with respect to the sky rather than 
with respect to the detector. There is only one star (SAO 17323, my = 8.8, rriK = 7.6) in the 
region of the QSO 1700 field which produces visible stray light artifacts in Channels 1 and 2 
(stray light paths are different for channels 3 and 4). However, the structure of these artifacts 
are different in Channels 1 and 2, and more important they appear in the lower coverage edges 
of the full images which are excluded from the power spectrum analysis. 

Interpretation of the signal witJi Limber equation 

More details related to this section can be found in ref. [2] of the main paper. When- 
ever CIB studies encompass relatively small parts of the sky (angular scales 9 < 1 radian) 
one can use Cartesian formulation of the Fourier analysis. If the fluctuation field, 5F{x), 
is a random variable, then it can be described by the moments of its probability distribu- 
tion function. The first non-trivial moment is the projected 2-dimensional correlation function 
C{9) = {5F{x + 9)5F{x)). The 2-dimensional power spectrum is P2(?) = where the 

average is performed over all phases. The correlation function and the power spectrum are a 
pair of 2-dimensional integral transforms, P2iq) = 27r C{9)Jo{q9)9d9. In the limit of small 
angles, the projected 2-dimensional correlation function is related to the CIB flux production 
rate, dF/dz, and the evolving 3-D power spectrum of galaxy clustering, P3{k), via the Limber 
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equation given by: 

C{e) = jdz(^ J_J{r;z)du (SI-2) 

where = c^{dt/ dzfv^ + d\9'^ and{(r) = S P3 (/c) jo (A;r)/c^(iA; is the two point correlation 
function for 3 -dimensional clustering with the power spectrum Pz{k). In the power spectrum 
formulation this equation can be written as (see e.g. refs. [2,12,13,21] of the main paper for 
details): 

where dA{z) is the comoving angular diameter distance and the integration is over the epoch 
of the sources contributing to the CIB. Multiplying both sides of this equation by /2t^ leads, 
after simple algebra, to eqs. (1),(2) in the main text of the paper. The fluctuation of the CIB on 
angular scale 27r/g can then be expressed as: 



^f^^Fc,BA{qd-/{{z))) (SI-4) 
where the suitably averaged redshift {z) is given by: 

In other words, the fractional fluctuation on angular scale vr/g in the CIB is given by the average 
value of the r.m.s. fluctuation from spatial clustering over a cylinder of length cto and diameter 
~ q~^d'2^{{z)) at the effective redshift given above. 
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Figure 1: Spectrum of CIB fluctuations Top: Power spectrum of signal minus noise from the 
SI Fig. 2 averaged over wide bins to increase signal-to-noise. The errors are N^^^"^ correspond- 
ing to the cosmic variance, where Nq is number of Fourier elements at the given g-bin. Solid 
lines show the shot noise from remaining galaxies fainter than the limiting magnitude in Table 
1. Filled circles and the darkest shade error bars correspond to the QSO1700 data, open circles 
and intermediate shade error bars to the EGS data and triangles with the lightest shade error 
bars and lines correspond to HZF data. Bottom panels: fluctuations, [^-^^]^^^, vs 271 /q for 
the QSO1700 data. Dashed lines estimate the contribution from ordinary galaxy populations 
and ACDM density field with At=5Gyr: the top dashed line shows the upper limit which as- 
sumes that their (high-2;) clustering pattern remained identical to that today at z=0 and the other 
dashed lines correspond to (2;)=1,3,5 from top to bottom. Solid lines show shot noise from re- 
maining ordinary galaxies. Dotted and dash-dot-dotted lines show the estimated Galactic cirrus 
and zodiacal light contributions respectively assuming the power spectra to be P{q) oc g" with 
n = —2 typical of cloud distributions. The observed cirrus power spectrum is a little steeper 

with n 2.5 to -3 ElSSlEIlEOll in which case the lines wiU have a slope of (3 + n) /2=0.25 

to 0.5. 
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Figure 2: Colour properties of clipped maps Estimates of colour between channels 1 and 2. 
Circles correspond to ^i^//^, averaged over the bins centered at these angles, diamonds to 

{S2{x)5i{x + 9)) / {5i{x)6i{x + 9)) for A^^cut = 4 and triangles to the same quantity averaged 
over maps with A^cut = 4, 3, 2.5, 2. Squares correspond to /34i//342 averaged over maps with 
iVeut = 4,3,2.5,2. 
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Figure 3: Contribution to CIB flux from Spitzer IRAC galaxy counts at 3.6 /im and 4.5 
/im. The value of dF/dm is evaluated from counts data and is a rapidly decreasing function 
of m of the form dF/dM oc exp(— 6m) with h ^ 0.4 (the least squares fits are shown with 
solid lines). (At the other two IRAC channels the uncertainty in the fall-off is much greater, but 
we show below that definite conclusions can already be reached from the data at 3.6 and 4.5 
ixm). Assuming that the functional form of dF/dm does not begin to rise appreciably at still 
fainter magnitudes gives the CIB contribution from galaxies fainter than mo of F{m > mo) = 
h~^dF / dm\mQ- For the QSO 1700 field, we eliminate galaxies brighter than mvega — 22.5 at 
3.6 /xm and mvega — 20.5 at 4.5 /xm corresponding to AB magnitudes ~ 25 and 24 respectively. 
The annotation in each panel show numbers which correspond to the extrapolated total flux 
from galaxies fainter than the magnitude limits in Table 1 using the least-squares fits shown 
with solid lines. One can expect that galactic spectra at the appropriate range of wavelengths 
are at most as steep as that of Vega, a star with a Rayleigh- Jeans black body spectral fall-off at 
these wavelengths. If so, they would generally have K — > leaving galaxies with KZ22 
(e.g. |41|). Galaxies at K> 20 have median redshift ^1 (f42'|), so the above argument would 
place the galaxies remaining at 4.5 fim at z^l and those remaining at 3.6 fim still farther out 
to higher z. Comparison with the Lyman Break Galaxy candidates in the same area shows that 
a substantial fraction of these galaxies are at 2; ~ 2 — 4 with 2; ~ 3 being the median redshift 
lfT4ll . Similarly, a substantial part of these galaxies are confirmed to be star-forming systems at 
z ~ 2.3 (ESI). 
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Maps are available from 



_ http://haiti.gsfc.nasa.g0v/kashlinsky/LIBRA/NATURE/_ 



4_as flo-SI.eps, fIb-SI.eps, fIc-SI.eps, fId-SI.eps 



Figure SI- 1 : Point source subtracted and clipped ( A^cut = 4, A^mask = 3) weighted intensity images 
of the QSO 1700 field for channels 1 through 4 (top to bottom). These positive images are aligned in 
the horizontal direction. They are scaled from [-2.0,4.0], [-2.2,4.4], [-9.2,18.4], and [-6.0,12] nW m"^ 
sr^^ respectively. These ranges are proportional to the standard deviation for each image. The regions 
of clipped sources are indicated by black areas which are set to 0.0 for further analysis. 



Mops ore ovalioble from 



http;//haiti.gsfc.nasa.gov/kashlinsky/LIBRA/NATURE/_ 



.as f2a-Sl.eps, f2b-Sl.eps, f2c-Sl.eps, f2d-Sl.eps 



Figure SI-2: Same as Fig. ISTTI but for iVcut = 2. 
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Figure SI-3: Figure resolution reduced to conform to astro-ph reuqirements. Original figure avail- 
able from http://haiti.gsfc.nasa.gov/kashlinsky/LIBRA/NATURE/ Signal and noise from QSO 1700 
field (top), the HZF (bottom - red lines) and the EGS fields (bottom - blue lines). The power spectrum 
of the residual diffuse emissions (+noise) is plotted with its errors and the noise is shown with open 
diamonds. The errors correspond to ^/Nq, where Nq is the number of Fourier elements that went into 
evaluating each value of P2{q)- 
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Figure SI-4: Autocorrelation function for the QSO 1700 field for A^cut = 2.5, 3, 4. Thick hnes corre- 
spond to A^mask = 3; thin lines to A'^mask = 5. The fraction of clipped pixels in the QSO 1700 data at 
iVcut = 4 and A^mask = 3 is between 25% (4.5 (Um) and 17 % (5.8 urn). Dotted lines show \/\C{9)\ on 
scales where C{6) < 0. 
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